AAS 15-259 


POSITION ESTIMATION USING IMAGE DERIVATIVE 

Daniele Mortari* * Francesco de Dilectis; and Renato Zanetti 


This paper describes an image processing algorithm to process Moon and/or Earth 
images. The theory presented is based on the fact that Moon hard edge points are 
characterized by the highest values of the image derivative. Outliers are eliminated 
by two sequential filters. Moon center and radius are then estimated by nonlinear 
least-squares using circular sigmoid functions. The proposed image processing 
has been applied and validated using real and synthetic Moon images. 

INTRODUCTION 

This paper completes a study, previously performed [1], to process cropped Moon and/or Earth 
images. The theory presented, which has been applied to real and synthetic Moon images, is based 
on the fact that the edge points are characterized by the highest values of the image derivative. 



Figure 1. Image Processing General Flowchart 


The content of this paper follows the general flowchart given in Fig. 1 which summarizes the main 
processes. These steps are explained in details in the following sections. The proposed algorithm 
has the following steps: 
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• The geometrical illumination conditions and the expected size of the observed Moon are first 
computed. These data allow to understand if the Moon is new, crescent, gibbous, or full. 


• The image differentiation is obtained with a fourth order Richardson extrapolation. An ana- 
lytical surface example is provided quantifying the accuracy gain obtained using Richardson 
extrapolation. 

• The pixels of the gradient image thus obtained, are sorted in descending mode. This means 
that the first set of pixels have high probability of being edge points. Outliers are then elimi- 
nated by two sequential filters. The first filter validates an edge pixel based on eigenanalysis 
of a small gradient box around the pixel and by investigation of the same box in the original 
image. The second filter is an improved version of the RANSAC algorithm. 

• The edge point selected are then used to obtain a first estimate of the Moon center and ra- 
dius. Taubin best fitting [6] (for circles) and the Fitzgibbons algorithm (for ellipses) are used. 
For the current application, the eccentricity of the Earth is so small that the ellipse angular 
orientation is a parameter ill defined to estimate (using noisy edge points). In addition, the 
presence of the atmosphere and high-altitude clouds make the problem even more inaccurate 
for Earth images. Therefore the orientation of the Earth ellipse is assumed known from the 
attitude of the camera. 

• The final estimation of Moon center and radius is performed by non-linear least-squares using 
a circular sigmoid function to describe the Moon edge transition. 

ILLUMINATION AND IMAGE PROCESSING PARAMETERS 

The image processing software require the computation of some internal parameters. With ref- 
erence to Fig. 2, these parameters depend on the Moon size on the imager, p r (cannot be smaller 
than a given value) and on the illumination geometry, pi (to define the separation between new and 
crescent Moon and between gibbous and full Moon). In this section, based on the geometry (Moon 
and Sun positions, J2000) and on a rough estimate of Spacecraft position (J2000), the conditions 
when the Moon is geometrically considered new, full, gibbous, or crescent are here derived. 


. Pi . 

k > i 



Figure 2. Expected Moon radius ( p r ) and illumination parameter (pi ) 


The observed Moon radius (in pixel) in the imager, p r , as well as the illumination size along the 
axis of symmetry, pi , dictate the values of important internal driving parameters used in the image 
processing code. 
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To estimate the full or new Moon cases, we use the fol- 
lowing input data. Let r s be Sun-to-Moon rays direction 
(unit- vector), r Q be (estimated) Spacecraft-to-Moon vector, 
and R m the Moon radius. From these data we can compute 


Pr = 


f Rn 




V 0 1 


(pixels) 


( 1 ) 


the illumination angle is then computed as 


cos P = r s • r Q 


( 2 ) 


from which P is derived. The angle a m can be derived from 
the equation 

\r Q \ sin<a m = R m (3) 



Figure 3. Full/New Moon geometry 


Full or new Moon 

while the value of the critical distance, r cr , can be obtained by 

r cr sin P — R m 

Therefore, 


if IrJ < r c 


and 


P < n / 2 then pi = 2 p r (Full Moon) 

> tv / 2 then p L = 0 (New Moon) 


(4) 


(5) 


Gibbous Moon 

Gibbous Moon will occur when 

\r 0 \ > r cr and 




( 6 ) 


When these conditions are satisfied then (see Fig. 4) the pi 
illumination parameter can be computed using the following 
equations 
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Crescent Moon will occur when 

7 r 

\r 0 \ > r cr and •& > - . (11) 

When these conditions are satisfied then (see Fig. 5) the pi 
illumination parameter can be computed using the following 
equations 

x 2 — \r Q \ 2 + — 2 \r Q \ Rm sin?? x (12) 

and 
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Crescent Moon 
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GRADIENT-BASED IMAGE PROCESSING ALGORITHM 

The gradient-based image processing approach is based on the fact that the part of the image with 
strongest variation of gray-tones is the Moon observable hard edge. Pixels with highest derivatives 
are the most likely to belong to the hard edge. Some of the selected bright gradient pixels do not 
belong to the Moon hard edge (outliers). These outliers are eliminated prior to using a best fitting 
algorithm to determine the sought geometrical properties. 


Image differentiation 

The gradient of a discrete function (the imager) is computed by numerical differentiation. In 
our case, the gradient is computed using four point central difference with a single Richardson 
extrapolation. 

The row and column four point central partial differences computed at pixel [r, c ] are 


8 [I(r + 1 , c) - I{r - 1 , c)] - [. I(r + 2, c) - I{r - 2, c)] 

12 h 

8[/(r, c + 1) — 7(r,c — 1)] - [/(r, c + 2) — /(r, c - 2)] 


(16) 


These derivatives are accurate with order /i 4 , where h is the pixel dimension. The image gradient is 
then computed as 


g'(r, c ) = yjg?{r, cj + 5 ' 2 (r, c) 


(17) 
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A single Richardson extrapolation requires computing Eq. (17) one time with step 2 h (getting 
g' 2h ) and one time with step h (getting g' h ). Therefore, the true derivative can be written as 


f 9 | true — 92h C (2^) 

\ 9 | true — 


Equating these two equation we obtain 


9h ~ 92h = c ^(2 4 - !) ch ^ = Wh ~ 32h) 


(18) 


(19) 


and therefore 


9 I true ~ 9h + 2 4 — 1 ^ 


( 20 ) 


An example of application of a single Richardson extrapolation is done using the six-hump Camel 
back function 
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Figure 6. Six-hump Camel back function 


Fig. 7 shows the accuracy of a single Richardson extrapolation 

In this figure the errors with respect to the true analytic values are show for the four cases: 
• Top-left: 2-point central differentiation and no Richardson extrapolation; 
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• Top-right: 4-point central differentiation and no Richardson extrapolation; 

• Bottom-left: 2-point central differentiation and one Richardson extrapolation; 

• Bottom-right: 4-point central differentiation and one Richardson extrapolation; 


These results clearly show that: 1) the accuracy provided by 4-point central differentiation and 
no Richardson extrapolation is equivalent to that provided by 2-point central differentiation and one 
Richardson extrapolation and 2) the 4-point central differentiation and one Richardson extrapolation 
provides almost the machine error accuracy. For this reason the 4-point central differentiation and 
one Richardson extrapolation is selected to evaluate which are the pixels with greatest gradient. 

The image gradient pixels are sorted in descending mode. The pixels with highest gradient values 
(edge pixels) are, therefore, concentrated as the first sorted elements. Unfortunately, some high 
gradient pixels (outliers) also appear on the Moon body (because of the illumination and local 
topology) or external (e.g., some bright stars or visible planets). Therefore, the outliers must be 
eliminated. Two sequential approaches have been implemented in order to eliminate the outliers. 
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Box-based outliers Identification 


The first approach to remove outliers is based on local analysis around the pixels being considered 
as potential edge pixels (pixels with the highest gradient values). With reference to Fig. 8 , this 
approach is based on using a mask size (the figure shows a 7 x 7 mask) whose size depend on the 
expected Moon dimension in pixels, p r . 

In Fig. 8 the central pixel (in blue) is the one being considered as potential pixel on the Moon 
edge. Using all the pixels of the gradient mask as lumped masses, the eigenvalues and eigenvectors 
of the inertia matrix are computed. The eigenvalues ratio (min/max) must be lower than a given 
threshold (e.g., 0.3), the eigenvector associated to the minimum eigenvalue (red line) in inclined by 
an angle 7 with respect to the row axis direction. In the figure, the pixels marked as green should 
both have high derivative values. The pixels marked as red, on the other hand should have very 
different gray tones. 



RANSAC for circle and ellipse 

A standard algorithm to eliminate outliers is called RANSAC which is an abbreviation for “RAN- 
dom SAmple Consensus”. The RANSAC algorithm [2] is an algorithm for robust fitting of models 
in the presence of many data outliers. The algorithm is very simple and can be applied ti lines, cir- 
cle, ellipse, etc. Given a fitting problem with n data points and unknown m parameters x , estimate 
the parameters. For example, for lines mn — 2, circles m = 3, and for ellipse m = 5. Starting with 
7V m in = n, the classic RANSAC algorithm perform Attest times the following steps: 

1 . selects m data points at random; 

2 . estimates the m parameter, meaning to find the curve passing through these m points; 

3. counts how many data items are located at a distance greater than a minimum value d m i n from 
the curve. They are N^. 

4. if < Nm in then set A^nm = A/^. 
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Figure 9. RANSAC flowchart 


One of the limitations of the RANSAC algorithm is the random procedure in which the points are 
selected. The random selection approach may persistently select outliers. To improve RANSAC, a 
technique (originally developed for Pyramid Star-Identification [3]) to select always different triads 
with indices continuously changing, has been adopted. 

Best fitting for Circle and Ellipse 

Once the pixels in the image have been selected, the algorithm must determine the geometric 
characteristics of the observed body, that is, radius and center for the Moon, considered a sphere, 
and axes and center for the Earth, considered an ellipsoid of revolution. Feature recognition for 
conic curves is based essentially on best fitting, which are typically grouped in two categories: 

Geometric fit tries to minimize the geometric distance between the data points and the conic. In 
general this method is considered the most accurate, however it does not admit closed solution and 
requires the use of iterative algorithms. 

Algebraic fit tries to minimize the “algebraic distance”, i.e. the implicit conic equation F(x,y) = 
0. Many methods based on this approach have been developed and they differ in accuracy, stability 
and complexity. They generally give results less accurate than the geometric fit, but on the other 
hand have closed form solutions and therefore can be implemented in a much simpler way. 

Taubin for circles 

Best fitting algorithms based on geometric fit are iterative methods requiring an initial guess, and 
their reliability and accuracy is ultimately solely dependent in the iterative scheme implemented. 
Typically, the Levenberg-Marquardt scheme is considered the most accurate. However, given the 
need to keep the computational load on-board at a minimum, algorithms based on algebraic fit seem 
a better choice, even if their accuracy is slightly worse. Among these, the Taubin [6] method is used 
and described in further detail. 

Given the circle implicit equation as 

a ( x 2 -\- y 2 ) + b x + cy + d — a z + b x + cy + d — 0 (22) 
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Method 

Geometric (G) 
Algebraic (A) 

Initial guess 
(Y/N) 

Notes 

Levenberg-Marquardt 

G 

Y 

Iterative, considered the “best” 

Trust Region 

G 

Y 

Improvement on the L-M 

Chernov-Lesort 

G 

Y 

Guaranteed convergence 

Kasa 

A 

N 

Fastest. Biased towards small radii 

Pratt 

A 

N 

Improvement on Kasa 

Taubin 

A 

N 

Reduces to a 3 x 3 eigenvalue problem 

Kukush-Markovsky-van Huffel 

A 

N 

Works best with big dataset 


Table 1. Summary of Best fitting algorithms for circle 


the Taubin Best Fit minimizes the function: 


T ( h _ n Ei( Q2: » + bxj + cm + df 

T ,C ’ ^Y(4a 2 Zi + 4abxi + Aacyi + b 2 + c 2 ) 

which is equivalent to minimize 

Ft ( a,b , c,d) = + bxi + cyi + d ) 2 

i 


subject to the constraint 

4a 2 z + 4afe x + 4 ac y + b 2 + c 2 = 1 

where the bar sign indicates simple average of the variables. 


(23) 


(24) 


(25) 


i i \ L \ L \ L 

X = ~ V Xi, y = ~y2vi, and 2 = - VV x 2 + y 2 ) = - V Zi . (26) 
2 = 1 2=1 2=1 2=1 

One important property of this method is that the results are independent of the reference frame. By 
reformulating this problem in matrix form, it becomes apparent how it is equivalent to a generalized 
eigenvalue problem. Indeed, define 


A= < 

"a" 

b 

> , z = 

~Z l X X y\ 1 

, and M = — Z T Z = 

zz zx zy z 
zx xx xy x 


c 

(dj 


z n x n y n 1_ 

n 

zy xy yy y 
z x y 1 


then the best fitting problem can be rewritten as 


where 


A T M A — 0 subject to A T N A — 1 , 

4 z 2x 2 y 0 
2z 1 0 0 

2y 0 10' 

0 0 0 0 


(28) 


(29) 
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Therefore, the optimization problem is equivalent to the following generalized eigenvalue problem: 


M A = r] N A. 


( 30 ) 


Since for each solution, [A, 77] 

A T MA = rjA T NA = r] ( 31 ) 

then to minimize A T M A the smallest positive eigenvalue should be chosen. Because the matrix N 
has rank three, such eigenvalue can be found by explicitly solving the cubic characteristic equation 
without the necessity to implement SVD or other more sophisticated algorithms. 


Fitzgibbons for ellipse 

The problem of best fitting for an ellipse is sensibly more complicated and somewhat less studied 
than the circle. Firstly, compared to the three parameters necessary to describe the circle (x, y,R), a 
generic ellipse is defined by 5 parameters (x, y , a, 6 , 6 ). Moreover, especially for very small values 
of the eccentricity, the inclination is very hard to discern and the error associated with it is several 
orders of magnitude greater than for the other parameters. 

Because the Earth itself, when considered an ellipsoid of revolution, has very small eccentricity, 
its planar projections are also almost circular, with ratio between the axes varying between 1 and 
0.9966 depending on the latitude of the observer. We thus expect the highest errors in inclination 
estimation. 

While geometric fit methods exist for the ellipse, they tend to not be very accurate and efficient 
due to the complicated formulation required. An non-linear Least Square approach, using the Ja- 
cobian of an ellipse in canonical form has been attempted; however, this method has proven itself 
not accurate enough for the requirements, and extremely slow to converge, even when advanced 
algorithms, like Levenberg-Marquardt are used. On the other hand, several studies ([ 4 ], [ 5 ]) have 
explored and improved the algebraic fit approach, with promising results. In particular, we have 
selected the Improved Fitzgibbon method, which reduces to solving a 3 x 3 eigenvalue problem. 

To show this, consider the implicit equation of a conic: 

/(a, 6 , c, d, e, /) = ax 2 + bxy + cy 2 + dx + ey + / = x T a = 0 . ( 32 ) 

For this equation to represent an ellipse, it must be b 2 — 4 ac < 0 , but with proper rescaling of the 
coefficient, which does not change the final result, it can be rewritten as b 2 — 4 ac = 1 . Therefore 
the problem can be rewritten in matrix form as: 

min a T D T Da = mina T 5 a subject to a T Ca — 1. ( 33 ) 

a a 

where 

0 0 2 0 0 0 

0 -1 0 0 0 0 

2 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

Solving this optimization problem leads to Sa = A Ca subject to a T Ca = 1, and the solution is 
the eigenvector corresponding to the minimum positive eigenvalue. However, solving the problem 


x\ xiyi y\ x x IJi 1 

D = : : : : : : , S = D T D, and C = 

x n x n y n y n x n y n 1^ 
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in this form leads to numerical instability. Instead, it is possible to decouple the quadratic part of 
D from the linear part, to obtain a simplified formulation that, together with being stable, has the 
added advantage of being of dimension 3 and thus the subsequent eigenanalysis can be performed 
without using SYD or similar. Therefore, by defining 



x\ Xiyi y\ 


xi yi 1 

D = [Z?i, D 2 \, where D\ — 


, and D 2 = 



_xl x n y n y\_ 


Xn Un 1_ 


and where 



'Si 
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s 2 ' 

S'3_ 


f Si 

= d\d x 


'0 

0 

4 

s = 

, where 

{ S 2 

= d\d 2 

and Ci = 

0 

-1 0 



l ^3 

= d t 2 d 2 


2 

0 

oj 


The final system of equation assumes the form: 


Mai = Ci 1 (S\ - S 2 SlS 2 )ai = A m i n ai 



and 



-s T 3 s 2ai . 


(35) 


(36) 


(37) 


(38) 


This version of the Fitzgibbon algorithm has been implemented with an explicit cubic solver for the 
characteristic equation. 


NUMERICAL TESTS 

To support the capability of the developed image processing code, three numerical examples of 
partially observed Moon are here presented. All these tests have been performed with no infor- 
mation about time, camera data, and observer, Moon, and Sun positions. The image derivative is 
obtained using the 4-point central approach with no Richardson extrapolation. These information 
clearly help to define the best values of the internal parameters. Since this help is missing, the 
internal parameters have been adopted with same values in all three tests. These values are: 


1. Maximum number of highest derivatives selected: N = 1500; 

2. Box-based outliers elimination: box size = 21 x 21; 

3. Box-based outliers elimination: derivative ratio >0.6; 

4. Box-based outliers elimination: corner ratio < 0.3; 

5. Box-based outliers elimination: eigenvalue ratio <0.3; 

6. RANSAC: minimum distance d m i n = 4 pixels; 

7. RANSAC: number of tests = 100. 
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Example #1: Real image, four- times cropped 

The original image of example #1, shown in Fig. 10, was intentionally cropped from a real image 
(taken from Earth) to validate the capability of the software to face a multiple cropped Moon images. 
Four times cropped is the also the maximum number of cropped segments possible. The four figures 
of Fig. 10 show the four main phases of the image processing. Once the image differentiation has 
been computed, the first 1,490 points with highest derivatives have been selected. This number 
differs from N = 1500 because the image is small and the number of points selected is also a 
function of the total number of pixels (N = 1500 is also the maximum number). 


Original linage 



Figure 10. Example #1: Original image 


1490 highest derivatives 21x21 Box approach: Removed 375 outliers 



100 200 300 

Number of tests = 100 Rac |ius = 21 8.0337 pixel. Worst residual = 2.9321 pixel 

Figure 11. Example #1: Image processing phases 


The box-based outliers elimination approach discarded 375 points based on the results shown 
in Fig. 12 while RANSAC eliminated additional 5 points. The resulting points are shown in the 
bottom right figure of Fig. 1 1 using blue markers. 

The bottom right figure of Fig. 1 1 shows the results obtained by Taubin fitting approach to 
estimate the circle parameters (Moon center and radius). The maximum distance from the estimated 
circle (worst residual) of the selected points was 2.9 pixels. The Taubin estimation gives a Moon 
center located at row 185.56 and column 187.34 and radius 218.03. Using this estimate, the set of 
points shown by red and blue markers in Fig. 13 over of Moon edge have been selected. These 
points are then used by CSF-LS to perform the nonlinear least-square estimation of Moon center 
and radius using circular sigmoid function. The new estimate implies a variation of —0.7, —0.5, 
and 0.4 pixels for row and column center coordinates and for radius, respectively. 
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CSF-LS Center (red) inside the image. Radius: 217.6452 pixel (A R = 0.38853) 



Column center: 187.8695 pixel ( A c = -0.53055) 


Figure 12. Box-based outliers statistics 


Figure 13. Example #1: CSF-LS final estimation 
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Example #2: Synthetic image cropped 

In the second example the test image is a synthetically generated image. All the relevant data and 
results are shown in Figs. 14, 15, and 16. 


Original Image Convoluted Image 



100 200 300 400 100 200 300 400 


Figure 14. Image convolution (3x3 mask) 


1500 highest derivatives 



100 200 300 400 

RANSAC: Removed 229 outliers 



100 200 300 400 

Number of tests = 100 


21x21 Box approach: Removed 911 outliers 



100 200 300 400 

Row = 189.0829 Col = 437.5094 



Radius = 89.0182 pixel. Worst residual = 4.9575 pixel 


Figure 15. Example #2: Image processing phases 


CSF-LS Center (red) inside the image. Radius: 89.0247 pixel (A R = -0.0064852) 



Column center: 437.1483 pixel ( A c = 0.3611) 


Figure 16. Example #2: CSF-LS final estimation 
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Example #3: Real image with center far outside the imager 


1500 highest derivatives 21x21 Box approach: Removed 35 outliers 



Number of tests - 100 Radius = 612.8029 pixel. Worst residual = 2.995 pixel 

Figure 17. Example #3: Original image 

Figure 18. Example #3: Image processing phases 


ERROR ANALYSIS 

This section is dedicated to the error analysis associated to the following sources: 

1. Estimated position error. This analysis quantifies how the error in the measured Moon 
radius (in pixel) affects the Spacecraft-to-Moon distance computation (in km). 

2. Estimated attitude error. This analysis quantifies how the error in the camera attitude 
knowledge (angle between quaternion) affects the Spacecraft-to-Moon direction (unit- vector). 

3. Image processing errors. This Montecarlo analysis shows the distribution of the Moon ra- 
dius and center estimated by the image processing software when the image is perturbed by 
Gaussian noise. This is done in the three cases of crescent, gibbous, and full Moon illumina- 
tion. 

Position and attitude uncertainty propagation 

This section investigates the error in the Moon estimated radius because of the position esti- 
mation error. That is, how the Spacecraft position estimate affects the measurement of the Moon 
radius in the imager. It also investigates the error in the Moon center direction because of the atti- 
tude estimation error. That is, since the Spacecraft position estimate requires the evaluation of the 
Spacecraft-to-Moon vector in J2000 while it is measured in the camera reference frame, how the 
camera to J2000 change of coordinate affect the Spacecraft-to-Moon direction? 
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Radius: 605.5698 pixel (A R = 7.233) 



Column center: -386.4042 pixel ( A c = -6.6663) 


Figure 19. Example #3: CSF-LS final estimation 


Let us consider the Moon radius (in pixel) uncertainty first. Moon radius is directly related to the 
Spacecraft-to-Moon distance. The Spacecraft-to-Moon vector is r om = r m — r where r m and r 
are the Moon and Spacecraft position vectors in J2000, respectively. Vector r m is simply derived 
from time while the Spacecraft position vector, r, it can be considered a random vector variable with 
mean f and standard deviation a r , that is, r ~ N(f, o r u ), where u is a random unit-vector variable 
uniformly distributed in the unit-radius sphere. Neglecting the ephemeris errors, the uncertainty in 
the Spacecraft vector position coincides with the uncertainty in the Spacecraft-to-Moon distance, 
O' D — O r - 

The Moon radius in pixels in the imager is provided by 


r 


d yz? 2 - w m 


(39) 


where d is the pixel dimension (in the same unit of the focal length, /), R m — 1, 737.5 km is 
the Moon radius, and D is the Spacecraft-to-Moon distance. This equation allows us to derive the 
uncertainty in the Moon radius in the imager 


C T r 


dr = [DJirn 

dD O d ~ d{D 2 -E? m f / 2 


(40) 
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The inverse of Eq. (41) 




d(z? 2 -i ^) 3 / 2 

f D R m 


CT rp 


(41) 


provides an important information about the accuracy required by the radius estimation. This equa- 
tion tells you how the radius estimate accuracy is affecting the Spacecraft-to-Moon vector length, 
and, consequently, how this estimate is affecting the estimation of the Spacecraft position. 



Figure 20. Calculated Moon distance error as a function of radius determination error 


Equation (41) has been plotted in Fig. 20 as a function of the Moon radius error and for four 
values of the Moon distance, D — 50,000 km, D = 150,000 km, D = 250,000 km, and D — 
350,00 km. 


The Moon direction uncertainty (in pixel) is caused by the uncertainty in the attitude knowledge. 
Consider the attitude q ~ A/"(q, 2cr q u) with mean q and standard deviation 2a q . Now the displace- 


ment with respect to the camera optical axis is c 



tan d, where d is the angle between estimated 


Moon center and earner optical axis. Therefore, the radial uncertainty in the Moon center direction 
is provided by 




/ 

d 


1 

cos 2 d 


cr n 


q 


(42) 
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Image processing error 

Image processing error quantification requires performing statistics using a substantial set of 
images with known truth data. Since no such ideal database is currently available, the analysis 
of the resulting estimated parameters (Moon radius and center) dispersion is performed by adding 
Gaussian noise to the same image. The Moon illumination scenario plays here a key role since 
when the Moon is fully illuminated the number of data (pixels) used to perform the estimation will 
be twice as much as in the partial illumination case (under the same other conditions: distance to 
the Moon and selected camera). For this reason two scenarios have been considered: 

1. Partially illuminated Moon (results provided in Fig. 21) 

2. Full illuminated Moon (results provided in Fig. 22) 


06-Mar- 20 13 05:24:30 CDT 



focal = 100 mm 


^ 0.3 
*3 

^ 0.2 
& 

s 0.1 
Z o 


- 0.1 

- 0.2 



- 0.2 0 0.2 
Col center error (pixel) 


1000 tests 


80 



Row center error (pixel) 


<t = 10 gray tone 


CSF-LS with 2823 pixels 



Radius error (pixel) 


Iterations 


Col center error (pixel) 


Figure 21. Image processing error results: partial illumination case 

The results of the tests done using partial illumination Moon are shown in Fig. 21. The image 
selected is real and it was taken on March 6, 2013 at 05:24:30 using a focal length / = 100 mm. The 
statistics of the results obtained in 1,000 tests by adding Gaussian noise with zero-mean value and 
standard deviation a — 10 gray tone are given. The code estimates the Moon radius and center (row 
and column) by least-squares using circular sigmoid function (CSF-LS). The results show the three 
parameters estimated as unbiased and with a maximum error of the order of 0.2 pixels (with respect 
to the original picture). The histogram of the number of iterations required by the least-squares to 
converge shows an average of 8 iterations. The central-upper plot in Fig. 21 shows the ellipsoid 
distribution of the Moon center error due to the fact that the estimation uses just data on one side of 
the Moon, on the illuminated part. 
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06-Mar-2013 05:24:30 CDT 


focal = 100 mm 


80 


1000 tests 




Col center error (pixel) Row center error (pixel) 



a = 10 gray tone 


80 



Radius error (pixel) 


CSF-LS with 2823 pixels 



Iterations Col center error (pixel) 


Figure 22. Image processing error results: full illumination case 


The results of the tests done using full illuminated Moon are shown in Fig. 22. Also in this case 
the image is real and was taken in Houston on February 25, 2013 at 22:33:00 CDT using a focal 
length / = 105 mm. The statistics of the results obtained in 1,000 tests by adding Gaussian noise 
with zero-mean value and standard deviation a = 10 gray tone are given. The code estimates the 
Moon radius and center (row and column) by least-squares using circular sigmoid function (CSF- 
LS). The results show the three parameters estimated as unbiased and with a maximum error of the 
order of 0.01 pixels, one order magnitude more accurate than what obtained in the partial illuminated 
case. The histogram of the number of iterations required by the least-squares to converge shows an 
average of 3 iterations. The central-upper plot in Fig. 21 shows the unbiased Gaussian distribution 
of the Moon center error. These more accurate results are due to the fact that the estimation uses 
data all around the Moon edge. 


CONCLUSIONS 

This paper presents an image processing algorithm capable of determining the apparent location 
of the Earth or moon and its apparent diameter. This information can then be used for autonomous 
onboard trajectory determination in cislunar space. Tow method of removing potential outliers are 
introduced in order to increase the algorithm robustness. Preliminary analysis of the algorithm’s 
performance is studied via processing both real and synthetic lunar images. 
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